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Gravitational wave emission from extreme mass ratio binaries (EMRBs) should be de- 
tectable by the joint NASA-ESA LISA project, spurring interest in analytical and numerical 
I methods for investigating EMRBs. We describe a discontinuous Galerkin (dG) method for 

• solving the distributionally forced 1+1 wave equations which arise when modeling EMRBs 

04 ' via the perturbation theory of Schwarzschild blackholes. Despite the presence of jump dis- 

continuities in the relevant polar and axial gravitational "master functions" , our dG method 
^ ' achieves global spectral accuracy, provided that we know the instantaneous position, veloc- 

. ity, and acceleration of the small particle. Here these variables are known, since we assume 

that the particle follows a timelike geodesic of the Schwarzschild geometry. We document the 
results of several numerical experiments testing our method, and in our concluding section 
discuss the possible inclusion of gravitational self-force effects. 
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PACS numbers: 04.25.Dm (Numerical Relativity), 02.70.Hm (Spectral Methods), 02.70.Jn (Collo- 
cation methods); AMS numbers: 65M70 (Spectral, collocation and related methods), 83-08 (Rela- 
tivity and gravitational theory, Computational methods), 83C57 (General relativity, Black holes). 



CN ■ I. INTRODUCTION 

> 

00 . An extreme mass ratio binary (EMRB) is a system comprised of small mass-mp "particle" 

' (possibly a main sequence star, neutron star, or stellar mass blackliole) orbiting a large mass-M 

blackhole, where the mass ratio fi = nip/M ^ 1. EMRB systems are expected to emit gravitational 



radiation in a low frequency band (10 to 10 Hz), and therefore offer the promise of detection by 



the joint NASA-ESA LISA project [H, Q]- A standard method for studying some EMRBs uses the 
C5 ' perturbation theory of Schwarzschild blackholes in an approximation which treats the particle as 

point -like and responsible for generating small metric perturbations which radiate away to infinity. 
These perturbations influence the trajectory of particle, resulting in deviation from geodesic motion. 



^ • Nevertheless, as a first and useful approximation, one may compute the emitted gravitational 

^ radiation, assuming that the particle worldline is a timelike geodesic in the Schwarzschild spacetime. 

More sophisticated approaches have included the effect of radiation reaction on the particle, usually 
through incorporation of a self-force (see Ref. 0] for a review) or a suitable approximation thereof 
[1]. Save for comments in the conclusion, this paper ignores self-force, although our methods may 
prove useful when it is included. 

A chief goal of gravitational wave signal analysis is to determine the spacetime structure of 
a given binary system from its experimentally measured waveform [3, IB] , a goal likely facilitated 
by numerical simulation. For the scenario we consider, simulation of ERMBs entails numerical 
evolution of Schwarzschild perturbations. The theory of such perturbations is well studied, and 
starts with pioneering investigations by Regge and Wheeler jg] and by Zerilli 0]. A detailed 
presentation of the subject is beyond the scope of our paper, and we point the reader to Refs. @, @] 
for modern accounts of the subject which supply the necessary background and references for our 
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work. Nevertheless, in order to provide some context, we now give a brief overview. We consider 
a small perturbation Sg^y of the background Schwarzschild metric g^y given in (j4j), where the 
perturbation satisfies the linearized Einstein equation. In our scenario the stress-energy tensor 
T^u given in ()C14p corresponds to a material point particle, and is therefore a distribution. The 
metric perturbation 5g^u is clearly tensorial; nevertheless, the perturbations can be reconstructed 
from a collection of scalar master functions. Remarkably, the master functions are governed by 
forced scalar wave equations with the following form:^ 

- dj^,^ + dl^.m - Vl{r)^lm = fir) [GUt, r)5{r - rp{t)) + F^t, r)6'{r - rp{t))] . (1) 



The coordinates here are the areal radius r, the Regge- Wheeler tortoise coordinate [10[ x = 
r + 2Mlog(^r/M — 1), and the time-dependent radial location rp{t) of the particle. Furthermore, 
Vi^r) is a potential, and ^£mit,r) is one of the master functions. Since master functions arise 
in a (tensor) spherical harmonic decomposition of the perturbative equations for dg^u, they carry 
multipole indices {£,m), where i > 2, \m\ < I. The distributional inhomogeneity on the right- 
hand side of stems from T^^u, and it involves Dirac delta functions, as well as the ordinary 
functions Fi^it,'''), G£mit,r), and /(r) = 1 — 2M/r. Perturbations of the Schwarzschild metric are 
characterized as either polar or axial, and each case corresponds to particular functions. The polar 
case corresponds to the Zerilli potential 0] 



{nir + 3M)2 



3M\ 9M^ / M\ 



(2) 



where = {£ + 2){£ — l)/2, and the Zerilli-Moncrief master function ^f^^. Appendix O provides 
explicit formulas for F™{t,r) and Gf^{t,r). The axial case corresponds to the Regge-Wheeler 
potential 



vrir) 



r 



(3) 



and for this case we choose to work with the Cunningham-Price-Moncrief master function ^^pm 
0]. Appendix ICl also gives explicit formulas for FpJ^^{t,r) and G^™(t,r). 

A number of numerical methods for solving ([T]) as an initial boundary value problem, and 
therefore modeling EMRBs in the time-domain, have appeared in the literature. In particular, we 



note Lousto's fourth-order algorithm based on spacetime integration of ([T]) and careful Taylor 
series arguments, and Sopuerta and Laguna's adaptive finite-element approach Jung, Khanna, 
and Nagle have applied a spectral collocation method to the perturbation equations for head-on 
collisions, using spectral filtering to handle the delta function terms [l^]- Most recently, Canizares 
and Sopuerta have proposed a multidomain spectral collocation method, with the particle location 



chosen between spectral elements 13|]. Clearly, the key difficulty to overcome is the distributional 
forcing; however, the problem should be amenable to a high-order accurate method, since — apart 
from possible transients — the solutions we seek to compute are everywhere smooth, except for 
jump discontinuities at the particle location. As a suitable high-order scheme for solving ([T]), 
we propose a discontinuous Galerkin (dG) method, and our approach shares some similarities 



^ We could instead work with the equation 

-df^em + 9^*£m - 14 (r)*^™ = gem{t)5{r - rp{t)) + Tt,n{t)5'{r - r-p(t)), 

where Gtmit) and J-em{t) depend only on t, and not on r. The relationships between Qimit) and J-em{t) and our 
Gim.{t,r) and Fem{t,r) follows from comparison between the right-hand sides of Eq. (lAlf) and the last equation. 
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with Refs. [l^, ™ particular we also ensure that the particle always lies at the interface 
between domain intervals. DG methods are widely used for wave-dominated problems, such as 
electromagnetic scattering [IJ, [l5| , and in nonlinear fluid dynamics 15|] . Our work is one of the 
first applications of dG methods to the modeling of gravitational waves (see also [3] ) , and the first 
dG computation of gravitational metric perturbations driven by a point-particle. Improving upon 



low-order methods, our method achieves global spectral accuracy (see also Refs. jl2l. Il3l|). 

This paper is organized as follows. Section [TI] provides further background necessary to under- 
stand the physical model. In particular, this section discusses the particle motion, the resulting 
jump conditions in the master functions, and a coordinate transformation adapted to the particle 
history. This background allows us to rewrite ([1]) as a first-order system which features only un- 
differentiated delta-functions in the forcing. Section IIIII describes our dG scheme as applied to the 
first-order system obtained in the previous section. Here we focus not only on the local represen- 
tation of the solutions, but also on how delta function terms are incorporated into the numerical 
flux function. Section IIVI documents the results of several experiments testing our method, and 
it compares our results with others found in the existing literature. In the concluding section, 
we discuss, in particular, the possible inclusion of radiation reaction. Several appendices collect 
technical results not given in the main part of the paper. 



II. PRELIMINARIES 

Throughout, we use both d and subscript notation for partial differentiation. For example, 
and are the same. We use an over-dot to denote d/dt differentiation, and sometimes a prime 
for differentiation by argument. The labels (^,m,CPM/ZM) are suppressed throughout. 



A. Particle motion 

In standard coordinates, the Schwarzschild line-element reads 



ds' = -fdt' + r'dr^ + r'{d9^ + sin^ 9d(f)'), (4) 

where t labels the preferred static-time slices and, as mentioned, r is areal radius. Owing to the 
spherical symmetry of the line-element, we may assume, without loss of generality, that the particle 



trajectory {rp(t),6p{t), 4>p(t)) = {rp{t),7r/2, (pp{t)) lies in the equatorial plane [l7|. Introducing the 
eccentricity constant e, semi-latus rectum constant p, and the parameterization rp{t) = pM/{l + 
ecosx(i))) we obtain the particle trajectory (f p it), (t>p( t)) by integration of the following system 
which describes timelike geodesic motion: [9|, \Vn, UM, llSl] 



# ^ (p - 2 - 2ecosx)(l + ecosx)^ , . 

dt Afp3/2[(p-2)2-4e2]^/^ 

dx (p — 2 — 2ecosx)(l + ecosx)^ [p — 6 — 2ecosx] 
dt [{p - 2)2 - 4e2] 



(5b) 



We use xit) rather than rp{t), since the former increases monotonically through radial turning 
points. In our scenario, integration of the system ([5]) is independent of ([1]). Therefore, we may 
view the particle path, and so the right-hand side of ([I]), as predetermined. We shall be interested 
in the parameter restriction < e < 1, for which the motion occurs between two turning points 
and the orbit is bounded. The periastron and apastron occur respectively at pM/{l + e) and 
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pM/{\ — e), and for e = the orbit is circular. Measured in coordinate time an eccentric orbit 
executes a radial period in time given by [l^ 
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Tr = C dx{l + ecosx)' 
Jo 
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2(3 + ecosx) 
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1 



2(1 + ecosx) 



1 -1 



(6) 
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When e 7^ 0, we average physical quantities over 4 radial periods, as defined by Eq. (j66p . 
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B. Jump conditions 

The forcing ([T]) induces jump conditions on the master function. Derived in Appendix El these 
are the following: 

{f^{t)-r'^{t))[m=U{t)F{t,r,{t)) (7a) 

2r,m M + i^pit) - mgpit)) M + - ^li^)) 

= fpit)G{t,rpit)) - gpit)Fit,rp{t)) - f pit) Frit, r pit)), (7b) 
where the subscript r in Frit, r) denotes partial differentiation with respect to the second slot, and 

fpit) = firpit)), gpit) = f'irpit)) / dtfpit) (8) 
are shorthands. In ([7]) our notation for a time-dependent jump is, for example, 

(9) 

Defining the particle velocity as Vpit) = Xpit) = Vpit)/ fp{t), we see that ([7a|) has the form 

fpit)il-v'pit))[m=Fit,rpit)), (10) 

confirming that the jump [[^']] is well-defined for a subluminal particle speed, \vp\ < 1. Therefore, 
we may safely make the substitution 



m it) ^ lim [M/(t, rpit) + e) - ^it, rpit) - e)] , 



fpit)Fit,rpit)) 



n{t)-riit) 

in all formulas which follow. Differentiation of (jlip gives 

2fpit)rpit)Fit,rpit)) [fpH) - fpit)gpit)] 



(11) 



+ 



im-r'pit))' 

gpit)rpit)Fit, rpit)) + /p(i)i^t(t, Vpit)) + /p(t)rp(t)F,(t, rpit)) 



(12) 



Finally, we may express (I7b|) as 

[[^r]] = [ - 2fpit)dt [M - if pit) - fpit)gpit)) m 

+ fpit)Git, rpit)) - gpit)Fit, Vpit)) - /p(t)F.(t, Vpit))] / if^t) - r^t)), (13) 

with the understanding that here and dt [[^']] respectively stand for (jll|) and (jl2p . Again 

note that fpit) — Tpit) > 1 for a subluminal particle speed, whence the jumps 9t [[^]] and [[^'r]] 
given by Eqs. ()12p and (|13p are finite. The formulas 

[M = dt M - rpit) [[^r]] , = Mt) (14) 

prove useful later. 
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C. Coordinate transformation adapted to particle history 

Now assume that x G [a, b] specifies the computational domain and the time-dependent particle 
location Xp = Xp{t) obeys a < Xp{t) < b, Vi. We enact the coordinate transformation 

t = X (15) 



u^p-a {b- Xp){ip- a) - {xp- a){b- ip, 

a + 7 4-a + 7- —--r ^ 4-a ?-?p, (16) 

ip-a {ip - a){b - ip){b - a) 



with the understanding that Xp = Xp{\) is explicitly time-dependent. The transformation obeys 
the following criteria: (i) t and A label the same time slices; (ii) x(A, |p) = Xp(A), with ^p = constant 
and a < E,p <b] (iii) x{\, a) = a and x(A, b) = b. We further require (iv) that the transformation 
is invertible on [a, 6]. This will only hold provided the point 

. _ iCp + a)(gp - XpjX)) + (xp(A) - a)(fc - gp) 

{^p-Xp{X)) ^''^ 

lies outside of the interval [a,b]. This is not a restriction of our method per se, and a coordinate 
transformation satisfying conditions (i) through (iv) may always be found. We have chosen to work 
with this one only for its simplicity. 
Differentiations of (jl6p yield 

dx (e-a)(6-e)xi(A) 



dX {^p-a){b-Cp) 

dx ^ (2^ -Cp- a){^p - XpjX)) + {xpjX) -a){b- Q 
dC {^p-a){b-Cp) 

d'^x _ 2{^p - XpjX)) 



(18) 
(19) 
(20) 



de {Cp-a)ib-^p)' 

and these expressions appear in later formulas. 

Under the coordinate transformation, the line-element ^ acquires a shift vector, 

ds'^ = -N^dX^ + L^{di + f3^dXf + r^ide'^ + sin^ 6d^'^). (21) 
Here = J^^^, L = f^^^dx/d^, with / understood as /(r(x(A, ^))). The shift vector is 

, ^ dx/dX ^ iC - a){b - Qx'pjX) 

^ dx/dC i2C-^p-a){Cp-XpiX)) + {xp{X)-a){b-Cp)' ^ ^ 

and we will also need 

_ {Ae + B^ + C)x'piX) 



9^ [m - ep - «)(ep - ^^(a)) + {xp{x) -a){b- e^)] 



2 ' 



(23) 



where A = 2{xp{X)-^p), B = 2{ab+(,'^-{a+b)xp{X)), and C = {a^+b'^)xp{X)-a{b-^p)^-b{a^+(,^). 
The velocity variable v = L(3^/N = dx/dX obeys 

t;(A,a) = 0, v{X,ip) = x'^{X), v{X,b) = {). (24) 

Since |xp(A)| < 1, we have |z;(A, .^)| < 1 uniformly in ^, assuming appropriately chosen ^p, o, and 
b. The vector field d/dX is not the Killing direction, and it does not point orthogonal to the 
constant-A slices. To relate the d/dX direction to the unit-normaP u of the slicing, first consider 



Here we use coordinate-free abstract notation for the vector fields u, u, and n. 
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g^^ = -Ar2 + = -Ar2(i _ y'i) = -{^Nj^f^ where 7 = (1 - 1;^)-^/^ is the relativistic factor. 

Therefore, the vector field 

u = -iN-^d/d\ (25) 

is normahzed, and one of its integral curves is the particle history. From standard formulas 

N~^{d/d\ - I3^d/d£) = 7"^^ - vL~^d/d(, (26) 

so that u = ju + v^n. Here g{n,n) = L~^g^^ = 1, whence n = L~^d/dS^ is a normalized spacelike 
vector field. These formulas show that the spacetime dependent parameter v determines a local 
boost in the tangent space of each spacetime point in the coordinate domain. At = this boost 
relates the slice normal u to the particle direction. 

D. Wave equation as a first— order system in (A, ^) 

Retaining the same letter ^ to denote the wave field ^(A, r(a;(A, ^))) in the new coordinates, 
we introduce the gradient $ = . Notice that = fdr^ in the old system corresponds 
to {dx/dO~^<^ in the new system. The following first-order system in the (A, ^) coordinates 



corresponds to the original second-order wave equation ([T]): 

5a = - n (27a) 

dxU = P^d^U - idx/dO~^d^[idx/dO~^^] + V{r)^ + Ji<5(e - Q (27b) 

dx^ = 95 (/?«$) - d^U + J2S{^ - (27c) 
where Eq. (j27h ) defines the variable H. The A-dependent functions 

Ji = -/3« [[H]] + idx/dO~' [[$]] , J2 = [[$]] + [[H]] . (28) 
implement the jump conditions collected in Section IIIBl where in terms of (I14p 

[[n]] = - [[^t]] , [M = (dx/dO [[^.]] . (29) 



The jumps (I28p can be recovered by integrating (I27p against a test function over the region (^p — 
e, + e)) performing an integration by parts, and taking the e — > 0"*" limit. Smooth terms vanish 
in the limit. Thus, the system (f27ll . with this choice of Ji and J2, is the first order form of ([I]) by 
construction. 



III. DISCONTINUOUS GALERKIN METHOD 



Following Ref. [15l |. this section describes the nodal discontinuous Galerkin (dG) method used 
to numerically solve (f27l) . Ultimately, we adopt a method-of-lines strategy, and here describe the 
relevant semi-discrete scheme which arises upon spatial approximation of (j27p by the dG method. 
In all numerical experiments considered later, we have carried out the temporal integration with 
an explicit fourth-order Runge-Kutta method, either one of low-storage [20| or the classical one. 
DG methods incorporate and build upon finite-element (FE), finite-volume (FV), and spectral 
methods, and in this section the reader will recognize the features our dG approach shares with 
these more traditional methods. For example, on each subdomain our approach features a weak 
formulation of Legendre collocation, and our technique for coupling subdomains draws on FV 
methods. 
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A. Local approximation of the system 

Our computational domain 0, is the closed ^-interval [a,b]. We cover Q with K > 1 non- 
overlapping intervals D'^ = [a^, b^], where a = a^, b = b^ , and b^~^ = a'' for k = 2, ■ ■ ■ , K . We 
further assume that the particle location S,p = b^p = a'^p'^^ lies at the endpoint shared by D^*' and 
Qfcp+i^ with 1 < kp < K. On each interval D'^, we approximate each component of the system 
vector (^',n, $) by a local interpolating polynomial of degree A^. For example, 

N 

*^(A,e) = 5^^(A,^,^K,^(e) (30) 

approximates where is the jth Lagrange polynomial belonging to D'^, 

4(o=n|-%- (31) 

Evidently, the polynomial interpolates ^' at the ^j. To define the nodes consider the 
mapping from the unit interval [—1, 1] to D'^, 

^^•(n) =a'= + i(l + n)(6*^-a^), (32) 

and the Legendre-Gauss-Lobatto (LGL) nodes Uj. The uj are the roots of the equation 

(l-n2)P^(n) = 0, (33) 

where Pn{u) is the iVth degree Legendre polynomial, and the physical nodes are simply = ^^(uj). 
In vector notation the approximation (130p takes the form 



*^(A,0 = *^,(A)^£^(e), (34) 

in terms of the column vectors 

^UX)=[^{X,eo)r-- ,*(A,d)]^ £'=(0 = [4(e),--- J^NiOf- (35) 

We also need to approximate by polynomials various products, for example (dx/dS,)^. Such ap- 
proximations are achieved through pointwise representations such as 

N 

ix^<^)U\0 = Y.x^{x,c^Mx,^^)e';{0. (36) 

j=0 

Here, and in what follows, we use the shorthands x^ = dx/dS^ and x^^ = d'^x/d^'^. Our vector 
notation for this example wih be (x^$)^(A,e) = {x^^)'l{X)'^ (.^{Cj- 

On each interval D'^ and for each solution component, we define local residuals, 

{R^t = d^^i-{(5^^)i + ni (37a) 

{Rnt = dx^i - d^{(5^Iit + i^d^f^^t + d^ixf^t + i^fx^^^t - {V^t (37b) 

[R^t = dx^i - d^P^^t + d^^l (37c) 

measuring the extent to which our approximations satisfy the original system of PDE. We define 
these residuals on open intervals {a^^b^) C D'^, but have assumed that the particle location = 
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b^p = a^p~^^ lies at an endpoint. Therefore, in the residuals (j37p we have not yet included the 
(5-function contributions appearing in (|27p . 

To motivate our derivation of a numerical approximation to (j27p , we first consider the kth inner 
product 

{u,v)c, ^ / d^uiCMO (38) 

J a'' 

and the expressions {(^j , {R^)h)D'=^ (^j i (-^n)^)D'=) ^^'^ (^j ) (-^'i')h)D'=- Namely, the inner products 
between the residual components (j37p and the jth Lagrange polynomial on D'^. We call the 
requirement that aZ/ these inner products {ij, (-R4')^)d'^> (^ji (-^n)^,)D*i i^j^ (-R*)h)D'^ vanish Vj the 
/cth Galerkin conditions. For now, we focus on the $ equation as a representative example, but 
will later also consider the ^ and IT equations. 

Enforcement of the Galerkin conditions on each D^' will not recover a meaningful global solution, 
since they provide no mechanism for coupling of the individual local solutions on the different 
intervals. Notice that, upon integration by parts, we may express the inner product as follows: 

(^•,(i?*)^)D. = r 9A$^:(A,e)+^f(0(/?«^)'(A,o-^f(e)n^(A,e)] 



- [(/3«$)^(A,e)-n^(A,e)]^,n0[l- (39) 
Therefore, in lieu of ([39|) with (-R$)^)d'= = 0, we enforce the equation 

a" 

This equation features numerical fluxes, -P^{X,a^){^'f^)* + (U^)* and -/3^(A, 6^)(<I>^)* + (H^)*, 
rather than boundary fluxes, -(/3^<I>)^(A, a*^) + n^(A, o'') and -(/3^$)^(A, 6*^) +n^(A, 5^^), thereby 
coupling adjacent subdomains. The numerical fluxes are determined by (as yet not chosen) func- 
tions 

(n^)* = (n^)*(n+,$+,n-,c|.-), ($1)* = ($^)*(n+ $+,n^,ci.-), (4i) 

where, for example, 11^ is an interior boundary value [either n^(A, a'') or n^(A, b^)] of the approxi- 
mation defined on D*', and 11^ is an exterior boundary value [either n^~^(A, b'^"^) or n^^^(A, a'^"'"^)] 
of the approximation defined on either D'^'"^ or D^~*"^. The fluxes ()4ip could also depend on 
but we will not need this extra generality. 

Let us now write the A^-|-l equations (j40p in matrix form, and also collect the corresponding 
matrix forms associated with approximation of the ^ and n equations. To write down these matrix 
forms, we first introduce the fcth mass and stijfness matrices, 

Mtj = f d^iUO^jiO, 4- = /' d^^fiO^j'iO- (42) 

These matrices belong to D^, and the corresponding matrices belonging to the reference interval 
[—1,1] are 



Mij = J duei{u)ej{u), Sij = J dui^{u)i'j{u), (43) 
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where ij{u) is the jth Lagrange polynomial determined by the LGL nodes uj on [—1,1]. These 
matrices are related by M^^- = ^{b^ — a^)Mij and S^j = Sij, whence only the reference matrices 
require computation and storage. In (|40p we now expand all polynomial approximations in the 
Lagrange polynomial basis as in (j30p . thereby obtaining 



M^dy^^i + {S^fiP^^t - [S'^VK = [M(A,0($'r - (n^)l^'(O t ■ (44) 
Now, as described in [3], the spectral collocation derivative matrix is 



(45) 



also given by = (Af^) ^S^. By the symmetry of the mass matrix, {D^y = {S^y {M^) ^. We 
then introduce the similarity-transformed matrix 

= M^D^{M^)-^, (46) 

tailored to obey {D\^y = {M^)~''{S^y . Applying {M^)~'^ to both sides of ([HD yields an equation 
for d\^\. The result and the corresponding equations for d\^\ and d\YL^ (derived via calculations 
similar to those presented above) are the following: 

a,*^-(/3«*)t + n^ = (47a) 

= {M'r' [ (/3«(A,o(n^)* - ^^-'(A, mir) ^Ho] (47b) 
d,^', + {Di,f{p^^t - {Di,fni = {M'r' [ (/3«(A, i){^ir - {niy) i\e)\ [1 . (47c) 

All adjacent vectors in these expressions, e. g. (/3^$)^, {V^)^, and (cc^'^cc^^^)^, should be inter- 
preted as a single vector obtained via component -by-component products. 

B. Numerical Flux 

To define the vector (/n, f^)^ of physical fluxes, we write ([27b .c) as 

( ^ ) ~'~ ( /'^ ) = lower order terms. (48) 
This equation determines the physical and numerical fluxes as follows: 



(49) 



The combinations of (n|)* and (<1>^)* which appear in (jlTb .c) are precisely — (/n)* and — (/|)*, 
as must be the case since these terms have arisen through integration by parts. In this subsection 
we construct the required boundary expressions for (/n)* and (/^)*. Our numerical flux must 
be robust, ensure stability, and be capable of handling the analytic discontinuities at the particle 
location. Numerical experiments suggest that inclusion of a Dirac delta function renders inadequate 
otherwise suitable numerical fluxes, such as the central and Lax-Priedrichs fluxes. However, as we 



10 



will see, a suitably modified upwind numerical flux successfully handles the delta functions in the 
system (I27p . recovering optimal convergence. We begin by constructing the standard upwind flux 
corresponding to no particle, and then incorporate the particle's effect into the flux through the 
addition of an extra term. 

An upwind numerical flux passes information across an interface in the direction of propagation. 
To construct the upwind numerical fluxes, we first diagonalize the matrix appearing in (I49p as 
follows: 



1 -/3«; I 



X 




Application of on the system vector (H, of fundamental fields yields the system vector 
(n + <I>/xg,n — ^/x^)'^ of characteristic fields. For our problem, the first characteristic field 
n + ^/x^ propagates rightward with speed — + x'^^ relative to the d/dX time axis, while the 
second characteristic field 11 — propagates leftward with speed — — x^^. Respectively, the 
upwind fluxes at a left endpoint a'^ (k ^ kp + 1) and at a right endpoint b'^ (k kp) then take the 
following forms: 

ifhr\ _^fO ,0 ^\^_,(U~\^(-P^ + xfO\^,(U-\ ^^^^^ 



Eqs. (|51b .b) formalize the intuitive concept behind the upwind numerical flux. In these equations 
triple-product matrices operate on the interior and exterior solution. The first matrix operation 
transforms the fields to characteristic fields, the second projects out one of the characteristic fields, 
and the third transforms back to the fundamental fields. As a result, information from a right- 
moving field, say, influences the subdomain to the right, but not the subdomain to the left. 

To achieve succinct expressions for the upwind flux which hold at both left and right endpoints, 
at each interface we define the average of a numerical variable and its numerical jump as 

{{<!>}} = + [[d,]]^^^ = n+^+ + n-^-. (52) 

Here n denotes the local outward-pointing normal of a subdomain and can be ±1. The numerical 
jump here is not a predetermined analytical jump as defined in Eq. ([9]), and it has a different sign 
convention. These definitions yield the following concise formulas (valid at left or right endpoints): 

(/n)* = {{-/3^n, + xf<!>H}} + \ [[xflih - xfa^^h]]^^^ (53a) 
(/I)* = {{n;, - (3^^h]] + \ Ixf^h - x5/3«n,]]^^^. (53b) 

At all interior endpoints (a'^ for k ^ 1, fcp + 1, and for k ^ kp, K) we will use this numerical fiux 
which is determined by the local numerical solutions. We also use this upwind form at a physical 
boundary (that is, or 5^), but in this case a boundary condition supplies the exterior solution. 

Turning now to the endpoints a^^'^'^ = b^f corresponding to the particle location, we modify 
the standard upwind fiux (j53p following the generalized discontinuous Galerkin method outlined 
in 



211 ]. Consider a Dirac delta function located at the interface between elements D p and D p 



and the weak form of the resulting system (j27p . The relevant new terms to consider have the form 
/ dCJi,2d{^-ip)ifiO, I dCJi,26{^-ip)£f^\0- (54) 



11 



Upon evaluation, each if these terms appears similar in form to a boundary flux. The discontinuous 
Galerkin method method provides a self-consistent way to evaluate these integrals and then add 
the results to the numerical flux. We only require the usual selection property of the delta function 
when integrated over the union 

Qfepy Dfcp+i^ and we are free to choose how the individual integrals 
over D'^p and D*'*'"'"^ contribute to the total integral. In fact, the dynamics of (j27p suggest a 
preferred distributional splitting. To see why, consider the scalar advection equation {d\ + vd^)u = 
J(^, A)(5(^ ~Cp)i with V > 0. Since this equation corresponds to rightward propagation, the natural 
choice for the associated distributional splitting of the delta function term is 

/ dCJ6{C-Cpyf{0=0, I d^J5{i-ip¥f^\£,) = J{ip,t)6^,r (55) 
For this case, notice that the delta function only "sees" a single Lagrange polynomial, namely 

k +1 

^0 (C) oi^ the rightward interval. 

To enact an upwind splitting of the delta functions appearing the system (I27p . we simply use 
the matrix already defined in ()50p to isolate the two propagating characteristic modes of the 
system. Consistent with propagation of these modes, at the particle location we modify the fluxes 
given in Eqs. (iSTk .b). 



left, modified 



U$ ^ / right, modified 




(56a) 
(56b) 



The correctness of this prescription can be see as follows. Integration of the system (j27p over the 
union D'^p U D'"'*'''"^ followed by a subsequent integration by parts on each interval generates the 
following boundary terms at the particle location (and on the right-hand side of the equal sign): 



/n 



(A,a'=P+i) 



/n 



The two physical fluxes in this equation of course cancel each other out, leaving only the vector 
{Ji, J2)^ ■ Our modifications ()56b .b) of the numerical flux are tailored to mimic this result. While 
the difference of the left /right numerical fluxes at the particle location will not, in general, cancel 
each other out (due to numerical error), notice that by subtracting ([56b ) from ([56b ) we generate 
precisely the vector (Ji, J2)^- This argument can be made more rigorous through an analysis 
based on integrating the two local numerical solutions on D^p and \^^p^^ against the Lagrange 
polynomials (.^{Cj ^-^d (0- Finally, using the general expressions ([53b ..b). we may likewise 
succinctly express the modified numerical flux at the particle location as 



modified V(/l) 
where either /c = /cp or /c = /cp + 1 in this equation 



C. Initial data and boundary conditions 



The issues of initial data and boundary conditions are not part of the dG method per se, but we 
must nevertheless specify both to complete our numerical scheme. We adopt trivial (zero) initial 
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k 






Im/3?-^(500) 


I 




-02 







— o.zoy ioD44UZOilj— 


HQ 


n 
U 


3 


-5.49064917188E- 


-03 





4 


-3.62410271081E- 


-03 





5 


-2.32805739548E- 


-03 





6 


-1.42584745587E- 


-03 





7 


-8.04688157035E- 


-04 





8 


-3.83719341654E- 


-04 





9 


-2.99532499571E- 


-03 


1.73407822255E-I 



k 


Rc/35 . (500) 




Im/35 , (500) 


\ 


-L.iitJ 1 OtyWiJVJi/ ( J. J_J 


-02 


Q 


o 
z 


— O.ZOUZy4Diyzill;— 


-Uo 


U 


3 


-5.48806353366E- 


-03 





4 


-3.62239165593E- 


-03 





5 


-2.32695433490E- 


-03 





6 


-1.42517041551E- 


-03 





7 


-8.04304980721E- 


-04 





8 


-3.83535015275E- 


-04 





9 


-2.99383340672E- 


-03 


1.73321233868E- 





Re72RW(5oo) 






k 


Re7|j^(500) 




Im7f ,,(500) 


1 


-8.36957985819E- 


-09 





1 


-8.35513276685E- 


-09 





2 


-2.95922379193E- 


-07 





2 


-2.95425498144E- 


-07 





3 


-2.97720676842E- 


-06 





3 


-2.97239482588E- 


-06 





4 


-8.13540247121E- 


-06 





4 


-8.12342297064E- 


-06 





5 


-1.40566197350E- 


-06 





5 


-1.40379108037E- 


-06 





6 


-5.02202428400E- 


-08 





6 


-5.01539234399E- 


-08 





7 


-1.01094068265E- 


-09 





7 


-1.00959570760E- 


-09 





8 


-7.70486047714E- 


-12 





8 


-7.69439666825E- 


-12 





9 


-2.99056309897E- 


-03 


1.73610608573E-03 


9 


-2.98758843820E- 


-03 


1.73437449497E- 



TABLE I: Compressed kernels for ^ = 2, n/[2M) = 500, e = 10"^°. There are d = 10 poles and strengths, 
and complex conjugation of the ninth entries gives the tenth entries. Zeros correspond to outputs from the 
compression algorithm which are less than lO"'^" in absolute value. 



data, and avoid the issue of an impulsively started problem by smoothly "switching on" the source 
terms, as discussed in Appendix O At the boundaries we impose outgoing radiation boundary 
conditions. Both potentials ()2l3p behave differently in the ^ —oo and — > oo limits, whence we 
treat the cases ^ = a and ^ = 6 differently. Since l-2Mr-i = 2Mr-^exp(-r/(2M))exp(x/(2M)), 
in the x — oo, r 2M~^ limit both potentials are exponentially small. Therefore, with a being 
sufficiently negative, |y^^'^(r)| is zero to machine precision when r corresponds to ^ ~ a, and as 
an excellent approximation we may use the Sommerfeld boundary condition 

{dt"^ -d^'^){\,a) = ^ U{\,a) + ^{X,a)/x^{X,a) = 0. (59) 

In the x,r ^ oo limit, both the Zerilli and Regge-Wheeler potentials (j2l3|) behave like V^^''^ = 
£{£+ l)r~^ + 0(r"^). Therefore, were we to adopt a naive Sommerfeld condition at ^ = 6, the slow 
fall-off of the potential would corrupt the benefits of our high-order accurate method. Instead, we 
implement the radiation boundary condition described in j22l |. 

-n(A, b) + $(A, b)/x^{X, b) = ^ nf^'^iX - X', rfe)^(A', b)dX', (60) 

''"b Jo 

where rh = r{x{X,b)) = r{b) and O^^'^ is a time-domain boundary kernel. As indicated, this 
kernel is different for the Regge-Wheeler (here spin-2) and Zerilli cases, although we suppress this 
dependence wherever possible. 

We approximate the time-domain boundary kernel 17^ ~ as a sum of exponentials 

^e{t,rb) = l^^e,k[t,rh), ^i^k{t,n) = — exp I — 1. (61) 

k=i ^ ^ 
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Solution * 





FIG. 1: vE'-COMPONENT OF THE SOLUTION. The 11 and <I> components are qualitatively similar. 



The parameters 7^^fc(rfe/(2M)) and Pi^k{f'b/{2M)) determine the approximation Ei{t,rh), and they 
depend on the Regge-Wheeler or Zerilh case, the orbital index i, and the dimensionless boundary 
radius rb/{2M). The approximation is designed so that its Laplace transform agrees with the 
transform of fl£ to relative supremum error e along the axis of imaginary Laplace frequency, and the 
the parameters 7/.^ and f3i^k are the outputs from the Alpert-Greengard-Hagstrom compression 
algorithm [j^ . 23]. Theoretically, e is a long-time bound on the relative convolution error in the 
time domain, and it measures the accuracy of the boundary condition. Table H] collects the £ = 2 
kernels for = lOOOM and e = lO"^''. We evolve the constituent pieces of the approximate 
convolution via temporal integration of the ODE 



d_ 
dX 



Hf,fc(A - A', n)^{}^, b)dX' + Ee,k{0, rb)^iX, b), (62) 



carrying out the integration along side, and coupled with, the numerical evolution of the system 
(j27p . With this boundary condition, we are free to choose essentially any boundary = 6, so long 
as it lies to the right of the source. Our outer radiation boundary condition is especially useful 
when studying eccentric orbits, for which one must average quantities over many periods. 



IV. NUMERICAL EXPERIMENTS AND RESULTS 
A. Forced ordinary 1+1 wave equation 

We consider two scenarios involving an exact solution of the distributionally forced 1+1 wave 
equation with no potential. The first adopts exact initial data, and the second trivial data in 
parallel with the blackhole perturbation problem (a setting where trivial data is often chosen). 
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FIG. 2: Temporal convergence of the linearly moving particle experiment. Errors have been 
computed relative to a uniformly spaced x-giid and over all fields. The dotted line is a least-squares fit of 
the data points (the round circles). 

1. Wave equation with exact initial data 

For a fixed velocity v obeying \v\ < 1, we consider 

-df"^ + 9^^' = cost5{x - vt) + icostJ'(x - vt). (63) 

Closed-form exact solutions for the real and imaginary parts of ^ are given in Appendix [Bl and 
we will check the convergence of our numerically generated solution against the complexification of 
these exact solutions. After expressing (I63p as a first order system and adopting our dG scheme, 
we obtain the same equations as in ()47p . except now with a zero potential vector V . Our domain is 
comprised of two subdomains: to the left of the particle location Xp{t) = vt, and to the right 
of Xp{t). At Xp{t), the interface between and D^, we use Eq. ([55]) for the numerical fluxes (7^)* 
and (/I)*. At the physical boundary points we choose fluxes which enforce simple Sommerfeld 
boundary conditions, 

n(A,a) + $(A,a)/x5(A,a) = 0, U{X,b) - ^{\,b)/x^{X,b) = 0. (64) 

For this flrst experiment, we take initial data from the exact solution. 

Working with the global domain [a, b] = [—5, 5], we choose v = 0.4 and the final time tp = 3.0. 
For these choices the critical ^ value (|17p always lies outside of the global domain, although clearly 
the example becomes pathological for a final time tp near 12.5 (when the particle crosses the outer 
boundary). Fig. [T] shows the ^ component of the solution vector, and the $ and 11 components 
also feature moving discontinuities. Fig. [2] documents the accuracy after several evolutions, each 
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FIG. 3: Spectral convergence of the linearly moving particle experiment. Again, errors have 
been computed relative to a uniformly spaced x-grid and over all fields. 

with = 26 points, performed with decreasing temporal resolution in order to exhibit the fourth- 
order accuracy of the temporal Runge-Kutta integration. To compute errors, we have used the 
polynomial representations of the two local solutions, each computed with respect to the coordinates 
(A,^), to interpolate onto a uniformly spaced x-grid with 256 points where errors have been 
calculated. Fig. [3] demonstrates the spectral convergence of our method for this problem. Here N 
is the number of points on each of the two subdomains, and for each N we have chosen a At to 
ensure stability. 



In this scenario we again consider (|63p. although now choosing u = (a fixed particle analogous 
to a circular orbit) and trivial initial data ^' = 11 = $ = 0. As before, we will compute errors 
over all fields, after interpolation onto a reference grid, and against the exact solution. Before 
computing errors for the ^ variable when adopting trivial initial data, we have found it necessary 
to adjust the mean of our numerically generated ^ in order to ensure that it equals the mean of 
the exact ^ . For the problem (163p with Sommerfeld boundary conditions (|64p . the ^ component 
of the solution vector is only determined up to an additive constant. Due to the presence of the 
potential, no such ambiguity is associated with solving ([T]). 

Our first test involves the minimal two domain set up. Since the problem is now impulsively 
started, we smooth the source functions as described in Eq. ()Cip . choosing to = 0, r = 3, and 
5 = 10. For these choices, the source is smoothly "switched on" (to machine precision) and is 
fully on by t = 3. Resolution of the transition requires relatively many points, and we have chosen 

= 61 on each subdomain. For the final time tp = 10, we demonstrate temporal convergence 
in the left panel of Fig. [H We note that, as indicated in the figure, convergence is abruptly lost 
without the smoother. However, even without the smoother, by adopting multiple subdomains 



2. Wave equation with trivial initial data 
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FIG. 4: Temporal convergence with trivial initial data. The left panel compares the two-domain 
experiment run with and without the smoother, denoted by circles and crosses respectively. The right panel 
corresponds to multiple subdomains and no smoother. As described in the text, on each subdomain wc have 
fixed the mean of the numerical ^' to the exact value before computing errors. 



we also recover convergence to the exact solution (of course assuming tp > 5, so that the initial 
incorrect profiles can fully propagate off the domain). Indeed, the right panel of Fig. U] documents 
the results for the same problem, but now without smoothing and 20 subdomains, each with N = 7 
points. We explain this observation by noting that for A'^ = 1 our dG method formally becomes a 
FV method. Therefore, many low-order elements corresponds to a more dissipative numerical flux, 
and the extra dissipation smooths the oscillations stemming from our impulsively started problem. 



B. Blackhole perturbations 

1. Experiment: Zerilli equation with radiation boundary conditions 

This experiment involves the i = 2, m = 2 polar problem and a circular orbit with p = 7.9456, 
M = 1, and nip = 1. We choose trivial initial data at to = 0, with a smoother defined by r = 10 
and 6 = 10. Integrating to final time tp = 90, we first generate an accurate reference solution 
^ref on the domain [—100, 100], using 65+55 subdomains (65 to the left of the particle and 55 to 
the right) with = 37 nodal points on each. Here and below, we choose the time step At to 
ensure stability. At both endpoints x = ±100 we place Sommerfeld boundary conditions on ^'rcf, 
as physically no radiation reaches the endpoints by the final time. 

The experiment is to generate a second numerical solution ^ on the shorter domain [—50,6], 
where b = 30+2 log (15 — 1) ~ 35.2. We again evolve to final time tp = 90, now with the convolution 
radiation boundary condition (j60p placed at the outer endpoint x = b. The relevant Zerilli kernel 



17 



is defined in Table II of Ref. [24]. This kernel corresponds to ri,/{2M) = 15 and the tolerance 
e = 10^^^. At the inner endpoint x = —50 we again adopt a Sommerfeld boundary condition. 
For 30+15 subdomains with 33 points on each, the corresponding ^ is then compared against the 
reference solution ^'ref in the Loo norm. After interpolation onto a uniformly spaced grid with 853 
points, we have found that — ^'rcflloo — 8.2314 x 10"-*^^. 



2. Results: circular orbits 

This subsection compares our numerical results for circular orbits to those obtained by other 
authors. For brevity we restrict ourselves to £ = 2, but note that our method maintains its 
performance for higher i. We have considered higher i values, and a fuller compilation of our 
results will appear elsewhere. For our simulations, we have chosen^ M = 1 = rup, with .^max = 
^^max = 1000+2 log(500 — l) ~ 1012 and ^min = Xrain = —200 as the outer and inner boundaries. We 
have used 45+200 subdomains, each with = 21 points, and a smoother (jCip defined by r = 1000 
and 5 = 0.0002. For these choices, we have integrated to tp = 2500 with time step At = 0.005. 
With these parameters we compute waveforms with a relative error of better than 10~®. Radiation 
boundary conditions (|60p have been enforced through Table [H Other parameters or non-uniformly 
placed subdomains may prove advantageous, but we have not explored all possibilities. 

We first describe what we have measured. The luminosities of gravitational energy and angular 
momentum across an arbitrarily large spherical surface are determined from the master functions 
+ x,r) and "^f^iu + x,r). We view the retarded time u = t — x as fixed, but with r,x 
arbitrarily large. Note that x ~ r, as r ^ oo. In the r — > oo limit we have the energy and 
angular momentum luminosities across an infinite-radius spherical surface given by (see [1,1^ and 
references therein) 

^ = E E = + (65a) 

£>2 m=~e ^ ^ 

f _ f f _ im (^ + 2)! ^ Y,CPM,t,CPM I ,f,ZM,t,ZM^ /f-f-ix 

^ - 2-^ ^irn, J^im - ^ (£ - 2)! ^^"^ +^em^em)- (o5b) 



£>2 m=-£ 



Here the overbar and dot denote complex conjugation and time differentiation respectively. The 
individual multipole contributions {Eim and Lf rn) to the total energy and angular momentum 



luminousities decay exponentially with i [1^, l2J, |25|. A few simplifications concerning Eim and 



are worth noting. First, due to the fact that the particle moves in the equatorial plane, the 
following conditions hold: £ + m even =^ xJ/CPM _ q ^^^^ £ -\-m odd =^ ij/^M _ establish 
these conditions, note, for example, that when £ + m is even the axial source terms and 
^fm^ are identically zero. Second, from the behavior of the master functions under the mapping 
m -m, we have Ei^rn = Ei_m and L£_„^ = Li_m [li|. 

We stress that Eqs. ([65k .b) hold at infinity, and practically one must devise a way to extract the 
waveforms at infinity from the finite computation domain. This problem has been solved for waves 



on flat spacetime by Abrahams and Evans 26|, |27[ • We will either simply "read off' waveforms at 



'^max = 1000 or use Abrahams-Evans flatspace extraction. For £ = 2 the procedure is as follows. 
We record a master scalar ^ at the outer boundary x = 6 as a time series, and then integrate 



^ By dividing Eq. fl} by Trip we can solve for the per-particle-mass perturbation ^ jm^ (from the coding stand- 
point, this is equivalent to setting = 1). Physical waveforms and other quantities can then be recovered via 
multiplication by appropriate powers of rrip. 
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Energy luminosity {E2m + -£'2, 


-m)/mp 




m 


dG, read off 


dG, extract 


FE 


FR 


FD 


1 


8.17530620 x lO"'^ 


8.1633 X lO-'^ 


8.1662 X 10"^ 


8.1633 X lO-'^ 


8.1623 X 10"^ 


2 


1.70685914 X lO""* 


1.7062 X 10"-* 


1.7064 X 10""* 


1.7063 X 10-"' 


1.7051 X 10"* 



Angular momentum luminosity (L2m + L2,-m) /rrv^ 


m 


dG, read off 


dG, extract 


FE 


FR 


FD 


1 
2 


1.83102416 X 10-5 
3.82285415 x 10"^ 


1.8283 X 10-5 
3.8215 X 10-^ 


1.8289 X 10-5 
3.8219 X 10-=* 


1.8283 X 10-5 
3.8215 X 10-^ 


1.8270 X 10-5 
3.8164 X 10-3 



TABLE II: £ = 2 luminosities for a circular orbit with (p, e) = (7.9456,0). 



Total i = 2 energy luminosity m.p ^ 5Zm=~2 (-'^Zm) 


Orbit parameters 


dG, read off 


dG, extract 


FR 


e = 0.18891539, p = 7.50477840 


2.59367 X 10-* 


2.59296 X 10-* 


2.59296 X 10-* 


e = 0.76412402, p = 8.75456059 


1.57146 X 10"* 


1.57120 X 10-* 


1.57131 X 10-* 



Total 1 = 2 angular momentum luminosity mp ^ 5I]m— -2(^2m) 


Orbit parameters 


dG, read off 


dG, extract 


FR 


e = 0.18891539, p = 7.50477840 


4.91165 X 10-3 


4.91018 X 10-3 


4.91016 X 10-3 


e = 0.76412402, p = 8.75456059 


2.09297 X 10-3 


2.09220 X 10-3 


2.09221 X 10-3 



TABLE III: Total £ = 2 luminosities for eccentric orbits. 



^f{t, b) ~ f{t — b) + 3/(t — b)b~^ + 3f{t — b)b~^ as if it were exact, thereby recovering the profile 
f{t) and its derivatives. We perform a similar extraction on 11. The Abrahams-Evans procedure 
is not exact for the perturbation equations we consider. Nevertheless, upon substitution of the 
approximate expansion f{t — x) + 3f{t — x)x~^ + 3f{t — x)x~^ into one of the (homogeneous) 
master equations ([1]), we find a residual which is ©(r""^ log(^r/M)) . Table HH compares our dG, 
circular-orbit, and i = 2 energy and angular momentum luminosities to results obtained by other 
numerical methods described in the literature. Such a comparison is not straightforward as the 
finite-element (FE) results of Sopuerta and Laguna involved reading off the master functions at 
X = 2000, while the finite-difference (FD) results of Martel [3] involved read-off at x = 1500 (here 
we always assume M = 1). The frequency-domain (FR) results of Poisson, as reported in [l8| . 
for the wave forms at infinity rely on the appropriate boundary value problems in the frequency 
domain, and of the three should afford the most direct comparisons. 



3. Results: eccentric orbits 

This subsection compares our numerical results for eccentric orbits to the frequency (FR) domain 
results of Tanaka et al [2^ (rather than Poisson's frequency domain results). We again choose 
45+200 subdomains, each with = 21 points, and At ~ 0.01. Due to the incommensurate radial 
Tr and azimuthal periods, we encounter the standard difficulty in obtaining measurements from 
eccentric-orbit simulations. Ideally, we would average measured luminosities over an infinite time, 
but will content ourselves with averaging over 4 radial cycles. Given a time series A{t), we compute 
its corresponding average as 

{A) ^ r dtAit), T2-T,= 4r,. (66) 

-t 2 — -t 1 JTi 
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FIG. 5: Orbital paths. The left panel shows one orbital period for {e,p) = (0.18891539,7.50477840). 
The right panel shows two orbital periods for {e,p) — (0.76412402, 8.75456059). In each case the dark inner 
circle is the horizon. Wc have used the (r, 0) system to construct these polar plots. 



Table Hill compares our total i = 2 angular momentum and energy luminosities to the frequency 



(FR) domain results of Ref. 28(] . In that reference the authors claim a relative numerical error of 



better than 10~^, which we have confirmed. We have retained enough significant digits in (e,p) to 
match the parameters (Ep, Lp) chosen in that reference. While we achieve relative errors of better 
than 10~^ for our averaged and extracted luminosities, we achieve single precision accuracy for our 
waveforms as a time series at x = b. Figure [5] exhibits the orbital paths for the two cases considered 
in this subsection, and Fig. [6] shows the corresponding waveforms. 



V. CONCLUSION 

We have presented a high-order accurate discontinuous Galerkin method for computing gravita- 
tional waveforms from extreme mass ratio binaries. Time-domain approaches for computing such 
waveforms have been hampered by the presence of distributional source terms (which include both 
a moving Dirac delta function and its derivative) in the governing master equations. By writing a 
master equation as a first order system, we have treated the source term physically through an ap- 
propriate modification to the numerical flux function. Our method maintains spectral convergence 
without requiring additional procedures (e.g. filtering), even pointwise in the immediate vicinity of 
the moving discontinuity. Through the use of convolution radiation boundary conditions, we have 
read-off waveforms at outer boundaries, thereby reducing computational cost without spoiling the 
high-order accuracy of our method. Accurate (read-off) waveforms, often with a relative error of 
better than 10~^, have been routinely observed in the course of our simulations. 

This work has assumed that the particle trajectory is a timelike geodesic of the Schwarzschild 
geometry. However, the gravitational perturbations induced by the particle will, in turn, affect 
the particle trajectory. Several existing techniques capture this effect, thereby incorporating more 
realistic inspiral (and possibly plunge) into the model. These include gravitational self-force con- 



structions 
calculations 
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jgl . 3ol . 31 1 (a representative, but far from exhaustive list), as well as post-Newtonian 



32l | and adiabatic approximations 
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FIG. 6: Waveforms for £ = 2, m = 2. The top panel shows the (e,p) = (0.76412402,8.75456059) 
extracted waveform, and the bottom panel the {e,p) = (0.18891539, 7.50477840) extracted waveform. Solid 
blue hnes and dashed red lines respectively correspond to real and imaginary parts. 



We believe that the central ideas of our approach might apply to many of these more sophis- 
ticated models. Computation of gravitational self-force requires that metric perturbations are 
regularized at the location of the particle. Mode-sum regularization has been carried out in the 
Lorenz gauge jso! ]. in an approach where the metric perturbations are described by the full coupled 
system of 10 PDEs rather than the simpler master equation description. Although a dG approach 
might certainly be applied in this setting, our methods are directly applicable to approaches based 
on master equations. Retaining the master equation description, Detweiler has argued for regular- 
izing gauge-invariant quantities [3]| . He has identified the appropriate quantities for quasi-circular 
orbits, and results based on these variables agree with corresponding Lorenz-gauge computations 
[s^ ]. Our numerical method should prove ideal for self-force computations which require that the 
metric perturbations are well resolved near the particle. In particular, we hope to use our method 
in tandem with self-force corrections based on regularization of g aug e-invariant quantities, at least 
for quasi-circular orbits. The post-Newtonian approach of Ref. [84] is also based on master equa- 
tions, and we might also follow that reference in order to include nonconservative effects in our 
simulations. 

In adapting our dG method to include, say, self-force effects we would surely encounter new 
difficulties. For example, due to dissipation associated with a self-force, the particle will inspiral, 
plunge, and merge with the blackhole. To handle all dynamical phases with our method, we would 
likely need to regrid at some point during the evolution. However, we believe this issue could be 
dealt with straightforwardly, given the robustness of the general method. 

We conclude with remarks on the applicability of our dG method to perturbations of the Kerr 
metric. Here, we consider only the scenario of a particle following a timelike geodesic of the Kerr 
geometry, although one might further consider gravitational self-force for this scenario as well. 
Now the relevant wave equation, the forced Teukolsky equation, is inherently 2+1 dimensional in 
the time-domain. In this case we would need to ensure that the particle always lies on an edge 
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between adjacent subdomains (in this case triangles). Clearly, this is a geometrically different 
problem, but Fan et al [211] have also considered 2+1 problems, and one might pursue the Kerr 
problem along similar lines. 
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APPENDIX A: JUMP CONDITIONS 

The derivation of the jump conditions ([7|) goes as follows. Using the selection properties of 5'{u) 
as a distribution, we first rewrite ([T]) as 

-dl^> + dl^ - V{r)^> = Mt)F{t, rp{t))d'{r - r,{t)) 

+ [/p(t)G(t,rp(t))-ffp(t)F(t,rp(t))-/p(t)F,(t,rp(t))]5(r-rp(t)), (Al) 

using the shorthand notations from Next, with L for "left" and R for "right", we let 

^{t, r) = ^^(t, r)9{rp{t) - r) + ^'^(t, r)0(r - rp{t)), (A2) 

where the step function 9{u) obeys 6{u) = for n < 0, and 9{u) = 1 for u > 0. We view the 
functions as everywhere satisfying the homogeneous PDF 

- + dl^^'^ - V{r)^^'^ = 0, (A3) 

even across the particle location rp(t). 

To complete our derivation of d?]), we calculate the distributional derivatives of ^ as given 
in ()A2p . insert them into ()Aip . and then compare terms. Using d/dx = fd/dr, the identity 
0'{u) = S{u), and the argument symmetry of S{u), we first compute 

= ^^e{rp{t) - r) + ^^e{r - rp{t)) - f^^^dir - rp{t)) + /^«5(r - rp{t)), (A4) 

where on the right-hand side we have switched to subscript notation for partial derivatives. The 
second ^-derivative of (jA2|) is then 

dl^ = ^^eirpit) - r) + ^^J{r - rp{t)) - 2/2^^'5(r - rp{t)) + 2f^^5ir - rp{t)) 
- ff^H{r - rp{t)) + ff'^'^dir - rp{t)) - /^^^^(r - rp{t)) + /^^^^(r - rp{t)). (A5) 

Expressed compactly, the last formula is 

dl^ = ^^,Mrp{t) -r)+ ^tO{r - rp{t)) + f^{t) [[%]]6{r - rp{t)) 

- fp{t)gp{t) [[^>]]5{r - rp{t)) + fp{t) [[^>]]5' {r - rp{t)), (A6) 

where the definition ^ is here [[^']] (*) = ^^{t,rp{t)) - ^^{t,rp{t)). To reach Eq. (jMj) from the 
previous line, we have used the selection properties of 5'{u). Next, we similarly compute 

= ^^Mrpit) -r)+ " rp{t)) + 2rpVl/f 5(r - rp{t)) - 2rp^f ,5(r - rp{t)) 

+ rp^^(5(r - rp{t)) - rp^f^5{r - rp{t)) - r2^^(5'(r - rp{t)) + r^^^^S'ir - rp{t)). (A7) 
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The last formula may be written in the succinct form 

dj^ = ^]:te{rp{t) - r) + ^^^(r - r,{t)) - 2rp{t) [[^t]] 5{r - rp{t)) 

- r,{t) [[^]]5(r - r,{t)) - fl{t) [[^.]]<5(r - r,{t)) + rl{t) [[^]5'{r - r,{t)), (A8) 

again by using the properties of 6'{u). Finally, with St [[5*]] = [[^t]] + rp(i) [[^'r]] , we rewrite the 
last expression as 

a^M/ = ^kO{rpit) -r)+ - rp{t)) - 2r^{t) {dt[[^]) 5{r - r^{t)) 

- f,{t) [[^]] 5{r - r,{t)) + rlit) [[*.]]<5(r - r,{t)) + rl{t) [[^]] 5'{r - r,{t)). (A9) 

Substitution of (IA6P and (lAOh into (lAlh . along with the fact that and solve the homogeneous 
PDE (jMI), then yields Eqs. ([TaD and (l7bl) . 



APPENDIX B: EXACT SOLUTIONS TO THE FORCED 1+1 WAVE EQUATION 

This appendix presents exact solutions to the distributionally forced 1+1 wave equation. Pre- 
cisely, we consider the equation 

- d'j-^ + dl^ = G{t)5{x - vt) + F{t)5'{x - vt), (Bl) 

where either F{t) = 0, G(t) = cost or F{t) = cost, G{t) = 0. In analyzing both cases, we make use 
of the following distributional identities: 

du\u\ = sgnti, duSgnu = 26{u), (sgnn)^ = 1, (B2) 

with sgnn = u/\u\ the sign function. Throughout, the particle location Xp{t) = vt has linear time 
dependence, with corresponding speed \v\ < 1. 



1. Solution for F{t) = 0, G{t) = cost 

Following analysis similar to that presented in Section IIIBI [or by substituting the correspon- 
dences Xp{t) = rp{t), X = r, f{r) = 1, f'{r) = 0, and f{t) = into the general jumps ([7])], we find 
the jump relations 

[[^L=.* = 0' P^n.=.t = 7' cost, [[dtn.=.t = -Vcost, (B3) 

where 7 = (1 — the usual relativistic factor. The particular solution 

^(t,x) = -isini?, {} = -f'^{t - XV -\x - vt\) (B4) 

to Eq. ()Bip possess the jumps listed in ()B3p . Using the identities ()B2p . let us verify that ()B4p 
indeed solves (IBip for F[t) = and G{t) = cost. Straightforward computation of the first and 
second order t-derivatives yields 

dt^ = -\'y'^[l + vsgn{x-vt)]cos-d (B5) 

d'^^ = + vsgTi{x - vt)f sin^ + v'^-i'^5{x - vt)cos^. (B6) 

while for the x-derivatives we similarly find 

d^-^ = \-i'^[v + sgii{x - viy^cosd (B7) 

= \-i\v + sg-n.{x - vt)fsm^ + -i'^5{x - vt)cos'&. (B8) 

Forming —d^^ + = 5{x — vt)cosi!}, we then appeal to the selection property of the delta 
function in order to reach the desired result, —df"^ + d^"^ = 6{x — vt) cost. 
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2. Solution for F{t) = cost, G{t) = 

The jump relations for this case are 

[[*L=.*=7'cost, [[d,^^^^^ = 2vjUint, [[d,^^^^^ = -{l + v'hUint. (B9) 
Now a particular solution to (|Bip is 



^'(t, j;) = 57^[z; + sgn(a; - tit)] cosi9, = 7^(t - - |x - (BIO) 

To verify that (jBTO]) indeed solves ([BT]) for = cost and G{t) = 0, we first calculate 

dt"^ = -l-f'^[2v + (1 + ?;^)sgn(x - 't;t)]sin?? - vj'^6{x - vt)cos'& (Bll) 
= v^^'^5'{x — vt) cos-!? + 7^[2i) + + sgn(x — vt)]5{x — vt) sin-f? 
- ^7^ [Sv + + (Su^ + 1) sgn(x - vt)] cos ^9, (B12) 

and then likewise compute 

a^^- = ^^^(x - ut) cos 1? + i7^[l + + 2v sgn(x - vt)] sin (B13) 
(9^'!' = 7^5' (x — ft) cos ■!? + 7^[3'L' + sgn(x — vt)]5{x — vt) sin •& 

- ^7^ [3t; + + (37;^ + 1) sgn(x - vt)] cos i?. (B14) 

Combination of Eqs. (jBTT]) and ([BT3]) yields 

- d"^-^ + 9^^' = 5'{x - vt) COST? + 7^[t> + sgn(x - ■ut)](5(x - vt) sin-d. (B15) 
By the selection properties of S'{u), we have 

6'{x — vt) COST? = 5'{x — vt) cost — 7^[f + sgn(x — vt)]6{x — vt) sint. (B16) 

Substituting this result into ()B15|) . using the selection property of S{u), and realizing that 
6{u) sgn u = by symmetry, we arrive at the desired result, — + = 6'{x — vt) cos t. 

APPENDIX C: SOURCE TERMS 

This appendix lists the specific functions Fijn{t,r) and Gimit,r) for our polar and axial cases. 
However, one practical difference between the functions listed below, and what we have used in 
our numerical simulations is the following. In accordance with our unphysical choice of vanishing 
initial data, we "switch on" the source terms smoothly via the following prescription: 

^ f ^[erf(\/5(t-to-T/2) + l]F^„(t,r) for to < t < to + r ^^^^ 
' I Fim{t, r) for t > to + r, 

and the same for G^m(t, r). Typically, the initial time to = 0, and the timescale r is much shorter 
than the final time of the run. Choosing suitable r and S, we ensure that our start-up is smooth 
to machine precision, and thereby avoid the troublesome nature of an impulsively started problem. 
Note that this prescription does initially affect the form of dtF£m{t,r). 

To both express and derive concrete expressions for the source terms, we rely on standard results 
for particle motion in the Schwarzschild geometry , [13]. As an alternative set to {e,p) discussed 
m Section UlAl we may instead work with {Ep,Lp), the physical particle energy and angular 
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momentum (both per unit mass). These constants of the motion are related to our original set by 



Li 



{p - 2)2 - 4e2 



^2 = - (C2) 

p - a - ^ p{p - 3 - e^) 

Let z^{t) = (t(r), r(r), ^(r), i?i)(r)) be the parameterization of the particle's four-trajectory in 
terms of proper time r. As before, set rp{t) = r(T(t)) for the radial coordinate of the particle 
expressed in terms of coordinate time, with similar expressions for Op{t) and 4'p{t). We assume 
that the coordinate system has been selected to ensure equatorial motion, 6p{t) = 7r/2. Then the 
four-velocity = dz^ /dr has components 



u 



r\2 



Ep/fir), K) 



f{r)il+Ll/r% 



u 



0, 



(C3) 



where these expressions follow from standard conservation arguments jlQ. Il7l|. 

Throughout this appendix, we use nip for the particle mass in order to avoid confusion with the 
azimuthal spherical harmonic index m. 



1. Zerilli— Moncrief (polar) source term 

To define the Zerilli-Moncrief source term, we first introduce the polar spherical harmonics 



Y 



Im 



Y 



im 



^ -.a 1 ^ ah 



^afe = ^-.ah H ^ y 7a6, 



(C4) 



where y^™'(0, 0) are the ordinary scalar harmonics 34 1, ^ah is the metric of the unit-radius round 
sphere, and a colon indicates covariant differentiation compatible with 7^^. The Zerilli-Moncrief 
source term is specified by^ 

/(r)F™(t,r) = e,(r)F^'"(t) (C5a) 
/(r)G™(t,r) = a,{T)r^{t) + 6,(r)y;"^(t) + c,{r)Yl^(t) + d,{T)Zl^{t). (C5b) 

Here, for example, Y^"^{t) = Y^^{tt/2, (t)p{t)). Moreover, the coefficients in (|C5P are given by 



a^(r 

ce{r 
di{r 
ei{r 



Svrmp /2(r) \ 6MEp A^(r) 



(l + ne) rAj{r) 
IQ-Krup f'^{r) L. 



En 



P r 

u 



(l + ni) r^Ai{r) Ep 

Snmp /3(r) 
{l + ne)r^Ae{r) Ep 

{£-2)lf{r) Ll 



-327rmK 



Svrmp f 



1 



{l + rii) Ai{r) Ep 




3M Ll 
l + n£ 1 — ^ { ne + 3 



7M 

r 



(C6a) 
(C6b) 
(C6c) 
(C6d) 
(C6e) 



where = {£ + 2){£ — l)/2 = A£{r) — 3M/r, and u'^ is determined by (|C3p and the sign of rp{t). 
Due to the u^' factor, we may not, strictly speaking, interpret ^^(r) as solely a function of r, but 
f{r)u^/Ep could also be reinterpreted as fp{t) and paired with ^^^"^(t). 



* Factors of /(r) are included here in order to have direct comparison with the same coefhcients listed in Refs. [9l[l8|. 
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2. Cunningham Price— Moncrief (gLxial) source term 

With €ab the unit-sphere Levi-Civita tensor such that ee^ = —s'mO, the axial spherical har- 
monics are 

oim ~,bc, \rlm olm olm (r^7\ 

'-'a —7 ^abJ^.c ' '->ab " '-'(a:6)' V*-"'^ 

and we express the Cunningham-Price-Moncrief source term as 

f{r)FP^''it,r) = Ce{r)Sl^{t) (C8a) 
/(r)GF™(t,r) = A,(r)5^-(t) (C8b) 

As before, {t) indicates evaluation on {9,<j)) = {tt/2, (pp{t)), and the coefficients in the above ex- 
pressions are as follows: 



A£{r) = 327rmp- 
Bi{r) = 32'jrnip- 
Ci{r) = 327rmp- 



(£-2)! 


P{r) Lp 


(£ + 2)! 


r2 El 


(£-2)! 


fir) Ll 


(£ + 2)! 


r3 El 


(£-2)! 


f{r) Lp 


(^ + 2)! 


r El 



(C9a) 
(C9b) 



L 



+ ^ 1 . (C9c) 



As before, we may not truly interpret Bi[r) as a function of r, but nevertheless keep this convenient 
notation. We note that our A(_[r) does not agree with the corresponding factor ui{r) quoted in 
Ref. 0]; however, we find that «f(r) = Ai{r) — C[{r). Due to this discrepancy, we present our 
derivation of (IC9D. 



3. Derivation of the axial source term 

Our goal is to establish formulas (|C8IC9p for the Cunningham-Price-Moncrief source term, 

^CPM ^ G^CPM^^^ _ ^ fCPM5'(^ _ ^^(i)). (CIO) 

Here and in what follows, we suppress m) indices wherever possible. Our starting point is Martel 
and Poisson's expression for 5^^^^ in (t, r) coordinates, 



^CPM _ _ 2r 



m+2) 



r^dtP' +fdrP' + ^P'^ . (Cll) 



This result appears in Appendix C of their expanded version of l8|] , where iSodd ii^ that reference is 
our The vector = {P\P'') is given by Eq. (5.10) of |J, 



but here with our index conventions and harmonics. Note that corresponds to of 
The stress-energy tensor for a point particle is 

T^'^ = mp [ dri-gy^/^u^u^S^x - z{t)), (C13) 
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here with g = — r^sin^0 the determinant of the background Schwarzschild metric We now 
change coordinates dr = {dT/dt)dt, integrate over t, and use = dt/dr, thereby finding 



P t 9 ■ 



5{r-rj,{t))5{e-ep{t))5{ck-cl)p{t)). 



(C14) 



Combination of Eqs. ()C12p and ()C14p . with the assumption of equatorial motion, then gives 



= ... . S^{7:/2, (t>p{t))5{r - rp(t)), 



+ 1) u 



(C15) 



Because we have not integrated over r, the four-velocity components u^^ here may be viewed either 
as functions solely of r, or solely of t upon replacing r by rp{t). Either viewpoint will yield the 
same derivatives dsP^ insofar as integration against test functions is concerned, and we view the 
components as depending on r. The delta functions of course depend on both r and t. Having 
identified which terms depend on r and i, we then calculate 



drP^ 
dtP'' 



Idirmp 

m + i 



idru'^)S^d{r-rp{t)) 



IGvrm^ 



(C16) 



£{i + l) n* 



(d^S^) - dr 



u 



5{r - rp{t)) - —S^drS{r - rp{t)) . (C17) 



To reach the last equation, we have replaced 4>p{t) by u'^/u^, which is permissible due to the 
presence of the accompanying delta function. Moreover, we have also made the replacement 



rp{t)5'{r - rp{t)) - ^5' {r - rp{t)) +[^] 5{r- rp{t)) 



w 



(C18) 



where the prime denotes partial r-differentiation. Finally, substitution of the last two results into 
Eq. (jCTT]) . along with Eq. jGH]) and the identity 



u'- 



2 fir) 
El r 



1- — ) + 1-/(0 



+ 



4M/(r) 



(C19) 



yields the desired results (jCSICOp . We have also used S^^{tt/2, (/)) = 5^S'<^(7r/2, i?!)), that is ordinary 
partial differentiation suffices in the equatorial plane. 
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